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Abstract 

The temperature and density conditions in the magnetized photosphere and chromosphere 
of the Sun lead to a very small degree of atomic ionization. In addition, at particular height, 
the magnetic field may be strong enough to give rise to a cyclotron frequency larger than 
the collisional frequency for some species, while for others the opposite may happen. These 
circumstances influence the collective behavior of the particles and some of the hypotheses 
of magnetohydrodynamics may be relaxed, giving rise to non-ideal MHD effects. In this 
paper we discuss our recent developments in modeling non-ideal plasma effects derived from 
the presence of a large amount of neutrals in the solar photosphere and the chromosphere, 
as well as observational consequences of these effects. 


1 Introduction 

The plasma in the photosphere and chromosphere of the Sun is often treated in the approx¬ 
imation of magnetohydrodynamics (MHD). MHD modeling has been very successful both in 
constructing idealized models of basic processes, and in providing realistic models of complex 
solar phenomena such as magneto-convection, formation of magnetic structures, flux emer¬ 
gence, wave propagation, reconnection, etc. [2113 [201 EH]. While in ideal MHD the plasma is 
supposed to be completely ionized, the solar plasma is only weakly ionized, and the ioniza¬ 
tion degree can drop as low as 10“^ in the upper photosphere. Together with the decrease of 
the collisional coupling with height, this may lead to the break of magnetohydrodynamic as¬ 
sumptions (depending on temporal and spatial scales of the phenomena under consideration) 
and lead to a series of non-ideal plasma effects. Since the degree of magnetization (the ratio 
of cyclotron to collision frequency of the different species) depends on the magnetic field, it 
is expected that non-ideal effects may play role in strongly magnetized photospheric regions 
(sunspots, flux tubes) and all over the chromosphere. 
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In order to take into account the non-ideal effects derived from the presence of neutrals, 
a ffuid-like description has often been used in the photosphere and the chromosphere. This 
is possible because the collisions are still strong enough to partially couple different plasma 
components. Many of single-ffuid and multi-ffuid models include only hydrogen or hydrogen- 
helium plasma, and a simplified treatment of radiative energy exchange (or even without 
taking radiative transfer aspects into account), which greatly simplifies the complexity of the 
problem. Such an approach has been frequently used in studies of different types of waves in 
the solar atmosphere [nmasHiiiaiis] and studies of instabilities laigEZ]. Idealized models 
of reconnection have been constructed in single-ffuid and two-ffuid approaches [271 [Ml is]- 
More complex models include those of magnetic ffux emergence [U [26] and of chromospheric 
thermal structure 1211 El EH]. The presence of neutrals, via collisions, provides a source of 
additional dissipation and energy supply. A weakly ionized plasma is able to diffuse through 
the magnetic field, partially breaking the frozen-in condition. The inffuence of neutrals 
results in creation or removal of cut-off frequencies for waves and instabilities, additional 
heating by dissipation of waves and currents, anisotropic dissipation of currents perpendicular 
to magnetic field facilitating the generation of potential structures, and other phenomena. 
Therefore, the inffuence of neutrals on the dynamics and energetics of solar plasma is proved 
to be important. 

With a constant increase of resolution capabilities of solar observations, both from 
space and from the ground, it is essential that new realistic models are built including non¬ 
ideal effects derived from the presence of neutrals. In this paper, we describe the efforts 
undertaken by our group in the development of an analytical theory and a numerical code 
for such modeling. We attempt to provide a self-consistent description of a multi-component 
multi-species solar plasma and its interaction with the radiation field. Direct observational 
detection of ion-neutral effects has also been attempted and its first results are described. 


2 Equations for partially ionized solar plasma 

We consider a plasma composed of a mixture of atoms of different atomic species. These 
atoms can be excited to excitation levels and/or to different ionization stages. Without 
loss of generality we limit the ionization stages to only two, as the number of multiply 
ionized ions is small in the regions of interest of the solar atmosphere. The interaction with 
radiation is included via excitation/deexcitation and ionization/recombination processes and 
photons are treated as another type of particles interacting with the rest of the mixture. The 
detailed derivation of the macroscopic equations of continuity, momentum and energy for 
each ionization-excitation level of species composing solar plasma components is given in 122 |. 
After summing up the motion equations over excitation states, a two-ffuid system of equations 
is derived that describes the behavior of neutral and ionized species. The electron and ion 
contributions are added up together to get a single equation for the charged component. This 
way, the following system of equations is obtained: 
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In the above equations pn — T^anT^an and pc = '^cxi'^cxi + n^eTTie are the total mass 
density of neutrals and charges; Un and are Uc their center of mass velocities; J is the total 
current density; = 3]9^/2 + Xan and e^i — ^Peil^^ + Zla Xm are the internal energies of 
neutrals and ions, where Xan and Xai are the ionization energies of species a. The summation 
in a goes over N chemical components. The definition of neutral, p^, and electron-ion, p^g, 
pressure tensors and heat flux vectors, and is done using a common system of reference 
for neutral and charged velocities, see 123 The heat flow vectors are then corrected for the 
potential ionization-recombination energy flux, e.g. ^ + Xan'^an- la this expression, 

the drift velocities of a neutral/ion of species a (e.g., Wan) are taken with respect to the 
central of mass velocity of all neutrals/ions, e.g. Wan — Uan — For the closure, one needs 
to provide the for collisional mass, momentum and energy exchange, Sn^ Rn and Mn> The 
mass source term Sn is due to imbalance between the ionization/recombination processes 
and depends on the number densities and on the collisional rate coefficients. The term Mn 
accounts for the energy exchange of neutrals through the elastic collisions with charges. Due 
to their complexity in a general case, we do not provide specific expressions for these terms. 
The collisional momentum exchange can be simplified to give 123: 

N N N 

Rfi ~ ~ Pei'^n ~ ^ ^ ~ Pii'^n ~ ^ ^ ^ ^ ^icxnp C^) 

P=l a=lP=1 

where z/gn^ are the collisional frequencies between electrons and neutrals of species /3 and 
are those between ions and neutrals of different species. 

The energy exchange by photons in Eqs. [^and|^is split into the contributions due to 
neutrals and charges to the radiative energy flux as follows: 

Fr^FI^Fr (8) 


where 

roo r roo r 
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In these definitions we have separated the contributions in the absorption kjy and emission 
coefficients related to neutrals and ions. The intensities Ijy are obtained from the complete 
radiative transfer equation, with total coefficients jjy and kj^. 

^=ju- Kh ( 10 ) 

a s 


To close the system dlHOh, we use the generalized Ohm’s law providing the electric held, E: 


E* = [E + UcXB] = 


Pe 
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This equation is obtained assuming stationary currents and neglecting second-order terms. 
From left to right the right-hand side terms in the Ohm’s law are: Ohmic, Hall, battery, and 
ambipolar term. The Ohm equation has similar form as the one for hydrogen plasma derived 
in e.g., Zaqarashvili et al. m except that we use Uc instead of Ui and include more complete 
expressions for the coefficients that depend on the collisional frequencies. 

The two-fluid formulation §§ is applied when the difference in behavior between 
neutrals and ions is larger than between the neutrals/ions of different kind themselves. The 
forces acting on neutrals and charges are different since the magnetic Lorentz force only 
affects charges, and the neutrals feel its action exclusively via collisions. Thus, the two-fluid 
approach is based on a stronger coupling between charged particles than between charged and 
neutral particles mms]. According to estimations of the magnetization factor of the solar 
atmosphere, the two-fluid approach is best valid at heights above 1000 km. Otherwise, when 
the collisional coupling of the plasma is strong enough, it is more convenient to use a single¬ 
fluid quasi-MHD approach i.e, including resistive terms that are not taken into account by 
the ideal MHD approximation. The conservation equations for the multi-species solar plasma 
become in this case: 

^ + V{pu) = 0 ( 12 ) 

^^ + Vipu0u) = JxB + pg-Vp (13) 

V (^u{e+^pu^)+ f>u ++VFr = JE* + pug (14) 

where we use a similar notation as for the two-fluid formulation. This system of equations is 
closed by the generalized Ohm’s law formulated for the single-fluid case, that reads as: 
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where G and an are given by 




( 16 ) 


N N N 

C^n = ^ Pel^enp + EE Pi^iatlp (17) 

/3=1 a=l /3=1 


and = Pn!P, are neutral and ion fractions. 

After applying Faraday’s and Ampere’s laws (neglecting Maxwell’s displacement cur¬ 
rent) to the Ohm’s equation formulated above, the induction equation is obtained for the 
evolution of the magnetic field: 
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(18) 


This equation closes the system of single-fluid equations. 


3 Extended Mancha code 


The single-fluid equations p^HU] and are solved by the extended version of the code 
ManchaSD. This is a non-ideal magnetohydrodynamic code with hyper diffusion algorithms 
and Cartesian grid written in Fortran 90. It solves the non-linear equations for perturbation in 
2D, 2.5D and 3D, and requires an arbitrary magneto-hydrostatic equilibrium to be explicitly 
removed from the equations. Spatial discretization is based on a six-order centre-difference 
scheme. The numerical solution of the system is advanced in time using an explicit fourth- 
order Runge-Kutta scheme. The code is fully MPI-parallelized using the distributed memory 
concept, and allowing full arbitrary 3D domain decomposition. A description of an early 
purely MHD version of this code can be found in m- Compared to the version described in 
that paper, the new code includes the non-ideal terms derived from the generalized Ohm’s 
law (Eq. 15), and the module for the calculation of non-ideal equation of state and the RT 
module that solves the Radiative Transfer equation (Eq. [Tol ). 


Order of magnitude estimates of the importance of the different terms in the generated 
single-ffuid Ohm’s law (Eqjl^ done in m for various atmospheric models show that only 
ambipolar. Hall and (potentially) battery terms can reach signiflcant values in the photo¬ 
sphere and chromosphere of the Sun. Therefore, in the current version of the code we only 
retain these three terms. This feature is partially described in [SUES]. 

The native set of variables used in the code is {p, pu^ e, B}. The pressure, p, is assumed 
to be a scalar and is obtained from e through the implementation of realistic Equation of 
State (EOS) for the solar mixture that takes into account the effects on ionization. This is 
done by assuming a non-ideal gas in thermodynamical equilibrium in the convection zone. 
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Figure 1: From left to right: initial distribution of the magnetic field; ion fraction am- 
bipolar diffusion coefficient rjA — (in m^ s“^); and the quantity J^rjA in a flux 

tube model. 


and an ideal partially ionized gas of a mixture including molecules in the photosphere, and is 
implemented as interpolations of lookup tables. An ideal gas with parametrized thermody¬ 
namical quantities is used as the first approximation elsewhere. We compute additionally the 
tables for electron pressure, so that the ionization fraction is obtained self-consistently with 
the rest of thermodynamical variables. We assume instantaneous ionization balance. The 
EOS implemented in that way provides the required precision while enormously reducing the 
computing time. 


The newly implemented RT module solves the non-gray Radiative Transfer equation 


(Eq. 10) assuming Local Thermodynamic Equilibrium (LTE). The wavelength dependence 
of the emission and absorption coefficients jy and ky is discretized by the opacity binning 
method m- The angle discretization is set prior to the code execution. Both discretizations 
are performed by sets of routines accompanying the code. The formal solver used in the RT 
module is based on the short-characteristics method m- The module is fully parallelized 
using MPI. The LTE approximation limits the application of the RT module to simulations 
of the top of the convective zone and the photosphere. In the optically thin corona, the code 
uses Newtonian cooling approximation. 


The code has passed standard tests on numerical performance and speed, with scaling 
tests in up to several thousand cores, showing its reasonable efficiency. The code is mod¬ 
ular, includes conditional compilation and a modern input/output file format (HDF5). An 
extended 2/2.5D version of the code uses Adaptive Mesh Refinement (AMR) implemented 
via the PARAMESH library. 


3.1 Chromospheric heating due to partial ionization effects 

Dissipation of currents due to the action of the ambipolar term provides an efficient source of 
heat. We performed 2.5D simulation of the chromospheric heating in an idealized atmospheric 
model composed of a set of flux tubes. We used a 2nd-order thin magnetic flux tube model 
constructed after [32], see Figure The model represents a horizontally infinite series of flux 
tubes that merge at some height in the chromosphere. This magnetic field configuration is 
non force-free. In order to evaluate the efficiency of the heating mechanism by ambipolar 
diffusion, we intentionally selected the temperature stratification of the flux tubes to be given 
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Figure 2: Simulations of chromospheric heating by ambipolar diffusion. Left: simulation 
snapshot at 800 s. The background color is temperature, white lines are magnetic field lines, 
arrows show the velocity field. For better visual comparison the temperature structure from 0 
to 0.42 Mm is kept unchanged from its initial value. Right: temperature as a function of height 
at horizontal position 0.6 Mm. Different lines are separated 0.5 sec in time; progressively 
more red colors indicate larger times till 800 sec since the start of the simulation. 


by a model HOLMUL [17], monotonically decreasing with height dictated by the conditions 
of radiative equilibrium. The ionization fraction in this model does not exceed 10“^, even in 
the chromosphere. The ambipolar diffusion coefficient, defined as tja = is orders 

of magnitude larger than the Ohmic one, reaching values of 10^^ m^ s“^ at a height of 2 Mm. 
Such values of r]A imply important current dissipation on very short time scales. 


The initial model is in MHS equilibrium and does not evolve without external per¬ 
turbation. As a perturbation we use ambipolar term in the induction (18) and energy (14) 
equations. Thanks to the Joule heating in the energy equation, JE^ the magnetic energy 
is efficiently converted into thermal energy, producingheat. This heat is balanced by the 


radiative cooling term, VFr. The left panel of Figure shows the simulation snapshot 800 
s after the introduction of the perturbation. The right panel of this figure shows the height 
dependence of the temperature at a fixed horizontal position inside the ffux tube for different 
time moments. The most important heating is achieved at the upper part of the domain, 
close to the tube borders, since the term responsible for the heating is orders of magnitude 
larger at these locations (right panel of Fig. [^. 


At the first few seconds of simulation, the chromospheric temperature significantly 
increases reaching up to 8 kK at the upper part of the domain, i.e. about 4 kK above its 
initial value. The time scales associated to the Joule heating are very short, of the order of 
seconds, similar to those of radiative cooling. Due to the balance between these two processes, 
after the initial heating, we observe damped temperature oscillations while converging to 
some constant value about 6 kK, forming hot chromosphere. Thus, we conclude that the 
Joule heating and radiative cooling terms in the energy equation can balance each other and 
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Figure 3: Time evolution of density in the simulations of Rayleigh-Taylor instability including ambipolar 
term (top) and without this term (bottom). The size of each snapshot is 1x1.8 Mm, the elapsed time is given 
at the bottom of each panel. Green lines are the projections of the magnetic field lines into the perturbation 
plane, the velocity field is indicated by arrows. 


lead to the chromospheric temperature rise. In a more complex situation, magneto convection 
simulations performed in [28] show that the action of the ambipolar diffusion allows to remove 
artificially cool bubbles produced by adiabatic plasma expansion in the chromosphere and 
to provide a more realistic model of its temperature structure. Therefore it is important to 
include ambipolar diffusion in further models of the chromosphere. 

3.2 Rayleigh-Taylor instability in solar prominences 

Plasma of solar prominences is yet another environment where the manifestation of partial 
ionization effects can be important. Typical values of temperature, derived from the Doppler 
width of spectral lines observed in prominences, are of the order of 8-9 kK [23110]. The 
ionization fraction is less known, but the values provided in the literature vary in the range 
=0.2-0.7 |2S|. The prominence-corona transition region (PCTR) is thought to be subject 
to tangential instabilities, as the Kelvin-Helmholtz instability (KHI) or the Rayleigh-Taylor 
instability (RTI) jU [151 ESI ES] . The linear theory has demonstrated that both instabilities 
are affected by the presence of neutrals in the prominence plasma [ZlEllST]. 

We performed 2.5D simulations of the non-linear phase of the RTI at the PCTR con¬ 
sidering partial ionization of the prominence plasma in the single-ffuid formulation. As initial 
setup we used a purely hydrostatic stratification of pressure and density for a given constant 
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Figure 4: Spectrum of the relative power of the “ambipolar” to “mhd” simulations of RTI 
with magnetic field normal to the perturbation plain {9 = 90°, top) and skewed from the 
plain at 0 = 89°(middle) and 9 = 88°(bottom). In the last two cases, the vertical dotted lines 
mark the cut-off wavelength of the instability, Ac = 38 and 155 km. 


temperature in the prominence (5 kK) and coronal (400 kK) parts of the domain. A constant 
horizontal magnetic field of 10 G was set over the whole domain, and its orientation with 
respect to the perturbation plane varied from one simulation to another. For these param¬ 
eters, the neutral fraction in the prominence makes = 0.9 and the ambipolar diffusion 
coefficient is r]A = 2.3 x 10^ m^ s“^. The instability is initiated by a multi-mode perturbation 
of the interface position. For each orientation of the magnetic field, we compare simula¬ 
tions with and without ambipolar term, see [23]. Figure shows a time series of simulation 
snapshots illustrating the development of the instability for the case with ambipolar term 
on (top) and off (bottom). The flows start the same, but then develop differences on small 
scales, since the ambipolar term introduces additional diffusion acting perpendicular to the 
magnetic field. These slight differences in the particular form of the turbulent flows can be 
observed at the end of the simulation. The ambipolar term is larger in the regions with 
low density, i.e. regions with, generally, low momentum pu. Therefore, the impact of such 
regions into the overall flow dynamics is not large. A similar behavior is observed in the 
simulations of multi-fluid turbulence in molecular clouds However, statistically, the 

flows in both cases are different, as is shown in Figure]^ In agreement with the linear theory, 
non-linear simulations demonstrate that the introduction of the ambipolar diffusion removes 
the cut-off wavelength for the growth rate of the instability and allows the small scales to 
develop. As the non-linear development of the instability is such that small scales merge and 
give rise to larger scales na [iHi sn, the larger growth rate of small scales leads to larger 
extreme velocities in the ambipolar case [23] . We also find up to 30% larger temperatures at 
the PCTR in the simulations with ambipolar term, as a result of the Joule dissipation. A 
significant drift momentum, defined as pjj — \/\piPn){ui — ^n) is present at the PCTR [23] . 

4 Observational manifestations of ion-neutral effects 

Direct observational confirmation of the uncoupled behavior of neutral and ionized species 
in the solar atmosphere is still missing. Order of magnitude calculations show that the 
scales at which ion-neutral effects become important are rather small. At maximum, spatial 
scales reach few kilometers, and temporal scales reach fraction of seconds, depending on 
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Figure 5: Slit-jaw image in Ha of the prominence observed on 11th of September, 2012. Black 
nearly horizontal line shows the location of the slit in one of the scans. 
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Figure 6: Top: amplitudes (left), widths (middle) and displacements (right) of the Call 4582 
A line as functions of the position along the slit and time. Displacements are given in km s“^. 
Bottom: same for the He I 10830 A line. 


the height in the atmosphere and on the magnetic field strength m- Direct observation 
of such small scales is not possible with current observational facilities. Nevertheless, if the 
differences in the behavior of ions and neutrals are sufficiently large in a small area inside the 
resolution element, it may be possible that their traces remain in the integrated signal even 
at the typical resolution of observations. This gives a chance of their detection in a specially 
dedicated observational campaign. Below we describe the preliminary results from two such 
campaigns. In one case, we target a prominence plasma. In the second case, the target was 
the Evershed flow in sunspots. 
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Figure 7: Left: velocity measured from Call line (red) and from He I line (blue) as a function 
of time at a fixed location. Notice the oscillatory behavior and slightly larger amplitudes in 
Ca at the velocity extremes. Right: velocities in both lines along the slit for a fixed time. 
Notice the lag between both velocities at the locations with larger gradients. 


4.1 Prominences 

The observations were done on 11th of September, 2012, using the Vacuum Tower Telescope 
(VTT) at the Observatorio del Teide on Tenerife. We targeted a prominence located at 
the limb close to the active region NO A A 11564 that was undergoing evolution during the 
observations. A slit-jaw image of the observed prominence is shown in Figure The data 
represent a time series of 10-steps scans (0.14" per step) of a small part of the prominence. 
It took 1.5 sec per scan position, a fixed position was repeatedly measured every 15 sec. We 
observed Call 8542 A and Hel 10830 A spectral lines, by simultaneously using the Tenerife 
Infrared Polarimeter II (TIP [6]) and the spectrograph of the VTT. Our primary objective 
was fast scanning in order to follow the rapid changes in the dynamics of the prominence. 
Therefore, no polarimetry was performed. The data acquisition started at 08:49 UT for 
almost 30 min, with the tq parameter staying around 10 cm. Standard reduction of the data 
was performed including the correction for the differential refraction for spatial alignment of 
the spectra taken at the two different wavelengths. The detailed description of the data and 
their analysis is given in [8]. 

The scientific analysis included a gaussian fit to the observed emission spectral lines, 
from which we extracted their amplitudes, widths and displacements [8]. The results of this 
fit are given in Figure]^ The upper raw of this figure gives the results for the Call line 
while the bottom one is for the Hel line. Due to the low sensitivity of the CCD at the 
wavelength of Call line the data for this line are much noisier than for Hel. Therefore, there 
are locations where we could not perform any fit of the Call line, or the quality of the fit was 
not adequate. In most of the locations both lines show simple shapes that can be interpreted 
as due to a single-component atmosphere along the line of sight. However, at some points 
the profiles show a clear two-component structure due to the appearance of the surges in the 
line of sight. In the preliminary analysis detailed here, we did not give any special treatment 
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Table 1: Parameters of the spectral lines in sunspot observations. Atomic data are taken 
from M- Formation heights are calculated in the penumbra models C and D from [36]. 


Elm Ion 

MV) 

FPL (eV) 

logy/ 

yeff 

He (km) 

Hd (km) 

I 

Fell 

4576.3 

2.84 

-2.976 

1.18 

449 

293 

0.69 

Fei 

4574.2 

3.21 

-2.500 

1.90 

280 

276 

0.48 

Fell 

4656.9 

2.89 

-3.643 

1.67 

162 

171 

0.41 

Fei 

4657.6 

2.85 

-2.900 

1.40 

277 

274 

0.41 

Fei 

5196.1 

4.26 

-0.493 

1.09 

336 

350 

0.69 

Fell 

5197.6 

3.23 

-2.348 

0.67 

539 

328 

0.72 

Fei 

5197.9 

4.30 

-1.640 

0.36 

216 

194 

0.39 

Fei 

5198.7 

2.22 

-2.135 

1.50 

398 

494 

0.79 

Fell 

5234.6 

3.22 

-2.279 

0.87 

553 

331 

0.72 

Fei 

5236.2 

4.19 

-1.497 

0.39 

235 

244 

0.38 

Fell 

6516.1 

2.89 

-3.432 

1.07 

281 

250 

0.46 

Fei 

6518.4 

2.83 

-2.460 

1.15 

327 

351 

0.51 


to these profiles. Figure shows that, outside of the locations where the Call line could 
not be fitted, the distribution of the amplitudes and widths of both lines in time and along 
the slit is rather similar. To first order, we take this as a confirmation that both lines form 
at the same location in the prominence and therefore measure the dynamics of the same 
plasma. It can be observed at the right panels of Fig. |^that indeed, the velocities derived 
from both lines are very similar. A more detailed look at the velocities is given in Figure 
It shows two selected examples of the time dependence at a fixed location (left) and spatial 
dependence at a fixed time (right). It can be observed that, while the velocities in Call and 
He I are very similar, differences appear at the locations of extreme velocities (left panel) or at 
the locations of more rapid changes in space where the velocity gradients are stronger (right 
panel). Having in mind the preliminary nature of our analysis, it may indicate that velocities 
of neutral and ionised atoms in prominences shows differences, partially attributing this to 
ion-neutral effects. 


4.2 Sunspots 

Observations were taken at the VTT during the morning of 7th of October, 2012. As a target 
we used an isolated regular-shaped sunspot belonging to the NO A A 11582 group located at 
coordinates (914'', -216") off solar disc center. By means of the spectrograph of the VTT 
we scanned the sunspot using several carefully selected spectral intervals. Each interval 
was selected from the the list of unblended spectral lines with a clean continuum provided 
in m under a criterium of containing at least one Fel and one Fell spectral line close 
in wavelength to fit on the same CCD. We selected pairs of lines with as similar formation 
heights as possible. The lines are listed in Table I. The scan step was 0.18" and 250 positions 
were taken. Different spectral intervals were scanned successively starting from 9:26 UT and 
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Figure 8: Continuum image and velocities of the Fel 5198.7 A and Fell 5197.6 A lines. The 
velocities are obtained after A-meter technique at the position around the core of each line. 
The circles mark the distances in the penumbra where we performed the fit to the velocities. 
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Figure 9: A-meter velocities as a function of height at the half radius of the visible penumbra, 
measured at different spectral intervals. 


it took about 5 min to complete each scan. The seeing conditions were exceptional during 
the morning, with the tq parameter varying around 15-20 cm. 

The data were reduced using a standard procedure and velocities were extracted by 
means of A-meter technique [39] . No absolute wavelength calibration was possible due to the 
absence of telluric spectral lines in the observed spectral intervals. We performed wavelength 
calibration by fitting the average quiet Sun spectrum to the FTS atlas. The 2D images of A- 
meter velocities and continuum intensity were rotated and corrected for the limb perspective. 
An example of the continuum intensity and velocities of the Fel 5198.7 A and Fell 5197.6 A 
lines are given in Figure]^ Our objective was to measure the velocity and the direction of 
the Evershed flow in the penumbra for neutrals and ions. For that we followed the procedure 
described in [3| and fitted the LOS velocity by a function 'yLOs (^5 0) = Vr{r) sinO coscj) + 
Vz{r) COS0, where r represents the radial distance along penumbra, 9 is the heliocentric angle, 
(j) is the azimuthal angle around the spot center, and Vr and Vz are the vertical and radial 
components of the velocity in the local reference frame. The fit was performed for the 
velocities at six radial distances in the penumbra located at various ranges of r, see Figure 

The A-meter technique allows to obtain the velocity at several positions along the line 
profile that correspond to different heights in the solar atmosphere. In order to compare Fel 
and Fell velocities, we need these heights to be reliably calculated. The formation height of 
spectral lines can strongly vary at the locations of bright and dark penumbra filaments. We 
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calculated the formation heights of the observed spectral lines in two semi-empirical models 
of sunspot penumbra, corresponding to a dark filament (model C) and a bright filament 
(model D) as provided in [36] at the observed heliocentric angle oi 9 — 77° (Shchukina, 
private communication). The height where the optical depth is equal to one at a given 
wavelength along the line profile was taken to be its formation height. Each A-meter velocity 
was assigned a height calculated this way in the models C and D and then averaged for a 
given r and A-meter position. For a reference. Table 1 provides heights corresponding to the 
line minimum. 

Figurej^shows the amplitude of the Evershed flow velocity v = ^/{v‘^-\-v‘l) as a function 
of height in the model D for all five spectral intervals separately for the Fe I lines (red) and 
Fell lines (black). The results are qualitatively similar for the heights from model C with a 
different temperature structure. In all the cases the outward velocity decreases with height 
according to the well known behavior of the Evershed flow. The Fe I lines show systematically 
larger velocities than Fell ones at all spectral intervals. Figure]^ shows the velocities at the 
half radius of the visible penumbra. The difference between the Fel and Fell velocities is 
maximum at radial distances close to the umbra-penumbra boundary and becomes almost 
negligible at distances close to the visible end of the penumbra. While the results of our 
analysis are subject to uncertainties due to the height of formation of spectral lines, it is 
intriguing that all spectral intervals independently show similar behavior with slightly larger 
velocity of neutral iron at the same photospheric heights. More extensive analysis of these 
data will be presented in our future publications. 


5 Concluding remarks 

This paper summarizes our recent effort in the theoretical and observational investigation of 
the effects of partial ionization of the solar plasma into its dynamics and energy balance. We 
have developed a self-consistent mathematical formulation for the two-fluid and single-fluid 
description of the partially ionized multi-species solar plasma and its coupling to the radiation 
field. The single-fluid equations are solved by the extended version of the 3D non-ideal MHD 
code Mancha. First results of our simulations show the high importance of the ion-neutral 
effects (ambipolar diffusion) for the heating of the magnetized solar chromosphere. The 
presence of neutrals is also shown to modify the stability criterium of the Rayleigh-Taylor 
instability at the border of solar prominences. Observational detection of ion-neutral effects 
has been attempted. Our preliminary analysis indicates possible differences in the velocities 
of ions and neutrals in solar prominences and in sunspot penumbra. 
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